program noiseadder2

    use iasing_osse
    use random_numbers_mix
    use netcdf
    
    implicit none
    ! Executable arguments
    character(len=128) ::  arg(5)

    ! Block properties definition
    integer, parameter :: BLOCKSIZE = 512
    integer :: nblocks
    integer :: jb, ja
    integer :: jamin, jamax, jalen

    ! Netcdf variables
    integer :: ionc, ncid
    integer :: iatm, natm, ilast
    integer :: open_mode
    type(r4a_ids) :: db_ids
    ! Netcdf data arrays
    real(JPR4)  :: wia(NIASI), wng(NIANG)
    real(JPR4), dimension(:,:), allocatable :: ria, rng, ngb

    ! Random numbers variables
    type(randomNumberStream) :: randstream
    !real(JPR8)      :: zrand_ia(NIASI*BLOCKSIZE)
    real(JPR8)      :: zrand_ng(NIANG*BLOCKSIZE)
    integer(JPI4)   :: kseed1, kseed2, kseedt
    integer, dimension(8) :: t8, rate8, cmax8

    ! Units change
    logical         :: lchangeunits = .false.

    !! Histogram variables
    !integer, parameter      :: nlimits = 81
    !real(JPR4), parameter   :: xlimits(nlimits) = &
    !                            (/(-5.+(2*5./(nlimits-1))*jb,jb=0,nlimits-1)/)
    !integer(JPI4)           :: histo(0:nlimits)

    ! Library includes
    include "r4a_initializer.interface"
    include "rad2brt_f.interface"
    include "brt2rad_f.interface"
    include "checknc.interface"
    include "histogram_f.interface"

    print*

    ! Initialisation of variables and modes
    open_mode = NF90_WRITE
    kseed1 =     524287
    kseed2 = 2147483647

    !! Get Last simulation
    !open(1,file='last',action='read',status='old')
    !read(1,*) ilast
    !close(1)
    !if (ilast.eq.0) then
    !    print*, "bad last value. Stop"
    !    stop -1
    !else
    !    print*, "ilast: ", ilast 
    !    ilast = ilast + 1
    !    print*, "ilast: ", ilast 
    !endif
    !!stop

    call get_command_argument(1,arg(1))
    print*, trim(arg(1))
    if (len_trim(arg(1)).eq.0) stop -1

    ! Open netcdf databank
    ionc = nf90_open(arg(1), open_mode, ncid)
    if (ionc /= NF90_NOERR) then
        print*, "Netcdf error: ", nf90_strerror(ionc)
        print*
        stop
    endif

    ! Initialise databank
    db_ids = r4a_initializer(ncid)

    ! Inquire length of sims
    ionc = nf90_inquire_dimension(ncid, db_ids%atm_did, len=natm)

    ! Get wavenumber arrays
    ionc = nf90_get_var(db_ids%dim%gid, db_ids%dim%wia, wia)
    ionc = nf90_get_var(db_ids%dim%gid, db_ids%dim%wng, wng)



    ! Determine number of blocks
    nblocks = (natm-1) / BLOCKSIZE + 1
    !print*, "natm:      ", natm
    !print*, "nblocks:   ", nblocks
    !print*, "BLOCKSIZE: ", BLOCKSIZE
    !print*

    ! allocate arrays
    allocate(ria(NIASI,BLOCKSIZE))
    allocate(rng(NIANG,BLOCKSIZE))
    allocate(ngb(NIANG,BLOCKSIZE))

    ! Couter format
100 format("END of block ",I6," of ",I6,", at ",I4.4,2("-",I2.2), &
            "_",2(I2.2,":"),I2.2,".",I3.3)
    ! Initilize counters
    call date_and_time(values=t8)
    write(*,100) 0, nblocks, t8(1:3), t8(5:8)

    ! Process blocks
    do jb = 1, nblocks

        ! Configure block limits
        jamin = BLOCKSIZE*(jb-1) + 1
        jamax = min(BLOCKSIZE*jb, natm)
        jalen = jamax - jamin + 1
        !print*, jamin, jamax, jalen, natm

        ! Get radiances
        call checknc(nf90_get_var(db_ids%rbt%gid, db_ids%rbt%ria,   &
                    ria(:,1:jalen), start=(/1,jamin/),              &
                    count=(/NIASI,jalen/) ))
        call checknc(nf90_get_var(db_ids%rbt%gid, db_ids%rbt%rng,   &
                    rng(:,1:jalen), start=(/1,jamin/),              &
                    count=(/NIANG,jalen/) ))
        if (ria(NIASI,1).lt.200) lchangeunits = .true. 


        if (lchangeunits) then
            ! Converting radiance to Brightness temperatures
            !print '(A,ES15.6)', "ria(1,1): ", ria(1,1)
            !print '(A,ES15.6)', "rng(1,1): ", rng(1,1)
            !print '(A,ES15.6)', "ngb(1,1): ", ngb(1,1)
            do ja = 1, jalen
                ria(:,ja) = rad2brt_f(NIASI, wia, ria(:,ja))
                rng(:,ja) = rad2brt_f(NIANG, wng, rng(:,ja))
            enddo
            !print '(A,ES15.6)', "ria(1,1): ", ria(1,1)
            !print '(A,ES15.6)', "rng(1,1): ", rng(1,1)
            !print '(A,ES15.6)', "ngb(1,1): ", ngb(1,1)
        endif

        ! Get random numbers
        call initialize_random_numbers(kseed1, randstream)
        call gaussian_distribution(zrand_ng, randstream)
        !call gaussian_distribution(zrand_ia, randstream)
        !print*, zrand_ia(1:10)
        !print*, zrand_ng(1:10)
        ! Update random seeds
        kseedt = modulo(int(kseed1,JPI8) + int(kseed2,JPI8), &
            int(huge(kseedt), JPI8))
        kseed2 = kseed1
        kseed1 = kseedt
        ! Write histogram
        !if (jb.eq.1) then
        !    call  histogram_f(nlimits, xlimits, NIANG*jalen, real(zrand_ng,4), histo)
        !    write(11,*) ja, " -Inf " , histo(0)
        !    do ja = 1, nlimits-1
        !        write(11,*) ja, (xlimits(ja)+xlimits(ja+1))*.5, histo(ja)
        !    enddo
        !    write(11,*) ja, " +Inf " , histo(nlimits)
        !endif

        ! Apply noise
        ria(:,1:jalen) = ria(:,1:jalen) +                       &
            reshape(zrand_ng(1:NIASI*jalen), (/NIASI,jalen/)) * &
            spread(NO_IASI, ncopies=jalen, dim=2)
        ngb(:,1:jalen) = rng(:,1:jalen) +                       &
            reshape(zrand_ng(1:NIANG*jalen), (/NIANG,jalen/)) * &
            spread(NO_IANGB, ncopies=jalen, dim=2)
        rng(:,1:jalen) = rng(:,1:jalen) +                       &
            reshape(zrand_ng(1:NIANG*jalen), (/NIANG,jalen/)) * &
            spread(NO_IANGA, ncopies=jalen, dim=2)

        if (lchangeunits) then
            !print '(A,ES15.6)', "ria(1,1): ", ria(1,1)
            !print '(A,ES15.6)', "rng(1,1): ", rng(1,1)
            !print '(A,ES15.6)', "ngb(1,1): ", ngb(1,1)
            ! Converting brightness temperatures to radiance units
            do ja = 1, jalen
                ria(:,ja) = brt2rad_f(NIASI, wia, ria(:,ja))
                rng(:,ja) = brt2rad_f(NIANG, wng, rng(:,ja))
                ngb(:,ja) = brt2rad_f(NIANG, wng, ngb(:,ja))
            enddo
            !print '(A,ES15.6)', "ria(1,1): ", ria(1,1)
            print '(A,ES15.6)', "rng(1,1): ", rng(1,1)
            print '(A,ES15.6)', "ngb(1,1): ", ngb(1,1)
        endif

        if (open_mode .eq. NF90_WRITE) then
            call checknc(nf90_put_var(db_ids%rbt%gid, db_ids%rbt%ria,   &
                ria(:,1:jalen), start=(/1,jamin/),                       &
                count=(/NIASI,jalen/) ))
            call checknc(nf90_put_var(db_ids%rbt%gid, db_ids%rbt%rng,   &
                rng(:,1:jalen), start=(/1,jamin/),                       &
                count=(/NIANG,jalen/) ))
            call checknc(nf90_put_var(db_ids%rbt%gid, db_ids%rbt%ngb,   &
                ngb(:,1:jalen), start=(/1,jamin/),                       &
                count=(/NIANG,jalen/) ))
        endif



        call date_and_time(values=t8)
        write(*,100) jb, nblocks, t8(1:3), t8(5:8)

        !exit
    enddo

    ! Deallocate variables
    deallocate(ria)
    deallocate(rng)
    deallocate(ngb)






    !do iatm = ilast, natm

    !    ! Get radiances
    !    ionc = nf90_get_var(db_ids%rbt%gid, db_ids%rbt%ria, ria, &
    !        start=(/1,iatm/), count=(/NIASI,1/) )
    !    ! Cycle if simulation not run
    !    if (ria(1) <= 0) cycle

    !    ionc = nf90_get_var(db_ids%rbt%gid, db_ids%rbt%rng, rng, &
    !        start=(/1,iatm/), count=(/NIANG,1/) )
    !    
    !    ! Convert radiances to brightness temperature
    !    if (ria(1).gt.1) then
    !        ! Change units from mWatts to Watts
    !        ria = ria/1E3
    !        rng = rng/1E3
    !    endif
    !    tia = rad2brt_f(NIASI, wia, ria)
    !    tng = rad2brt_f(NIANG, wng, rng)

    !    !print*, " -- ria: ", ria(1)/1E3
    !    !print*, " -- tia: ", tia(1)

    !    ! Initialise random seeds
    !    kseedt = modulo(int(kseed1,JPI8) + int(kseed2,JPI8), &
    !        int(huge(kseedt), JPI8))
    !    kseed2 = kseed1
    !    kseed1 = kseedt
    !    call initialize_random_numbers(kseed1, randstream)

    !    ! Get random seeds and modify brightness temperatures
    !    call gaussian_distribution(zrand(1:NIASI), randstream)
    !    tia  = tia + real(zrand(1:NIASI),4)*NO_IASI
    !    call gaussian_distribution(zrand(1:NIANG), randstream)
    !    tnga = tng + real(zrand,4)*NO_IANGA
    !    call gaussian_distribution(zrand(1:NIANG), randstream)
    !    tngb = tng + real(zrand,4)*NO_IANGB
    !    
    !    ! Return to radiance space 
    !    !print*, " -- tia: ", tia(1)
    !    ria  = brt2rad_f(NIASI, wia, tia)
    !    !print*, " -- ria: ", ria(1)
    !    rnga = brt2rad_f(NIANG, wng, tnga)
    !    rngb = brt2rad_f(NIANG, wng, tngb)

    !    if (open_mode .eq. NF90_WRITE) then
    !        ionc = nf90_put_var(db_ids%rbt%gid, db_ids%rbt%ria, ria, &
    !            start=(/1,iatm/), count=(/NIASI,1/))
    !        ionc = nf90_put_var(db_ids%rbt%gid, db_ids%rbt%rng, rnga, &
    !            start=(/1,iatm/), count=(/NIANG,1/))
    !        ionc = nf90_put_var(db_ids%rbt%gid, db_ids%rbt%ngb, rngb, &
    !            start=(/1,iatm/), count=(/NIANG,1/))
    !    endif
    !    print*, iatm

    !    !if (iatm==100000) exit
    !enddo



end program noiseadder2
